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Abstract 

We study thermal transport in a chain of coupled atoms, which can vibrate in longitudinal as 
well as transverse directions. The particles interact through anharmonic potentials upto cubic 
order. The problem is treated quantum mechanically. We first calculate the phonon frequencies 
self-consistently taking into account the anharmonic interactions. We show that for all the modes, 
frequencies must have linear dispersion with wave-vector q for small q irrespective of their bare 
dispersions. We then calculate the phonon relaxation rates Ti(q), where i is the polarization index 
of the mode, in a self-consistent approximation based on second order perturbation diagrams. We 
find that the relaxation rate for the longitudinal phonon, T x (g) oc q 3 / 2 , while that for the transverse 
phonon T y (q) oc q 2 . The consequence of these results on the thermal conductivity k(N) of a chain 
of N particles is that k(N) oc N 1 / 2 . 

PACS numbers: 44.10.+i,05.70.Ln,05.60.Gg,66.25.+g 
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I. INTRODUCTION 



In low dimensional systems, the description of transport in terms of linear phenomeno- 
logical laws, that works well in three dimensions, may not be possible. One example is 
the Fourier law which relates the thermal current density j to the thermal gradient VT by 
j = —kVT, where the k is the thermal conductivity of the material. Theoretical investi- 
gations of a number of translation invariant one-dimensional systems show that the Fourier 
law is not followed in these systems The ratio of the current density to the thermal 
gradient is no longer independent of the system size. Nevertheless this ratio, again called 
'thermal conductivity', is still useful in characterizing anomalous transport. 

The calculation of 'thermal conductivity' in these systems can be problematic. Ideally 
one would like to calculate the non-equilibrium stationary state (NSS) measure of the system 
and calculate the expectation values of j and VT using the NSS measure. This requires 
the existence of local thermal equilibrium (LTE) so that a local temperature T is properly 
defined. A method generally used when considering a one dimensional system with its 
ends at two different temperatures T\ and T 2 , is to calculate directly the current j and the 
quantity k(L) = jL/ (T 2 —T{), L being the length of the system. Another approach is to use 
the Kubo-Green formula which involves the time integral of the equilibrium current- current 
correlation function. The first formula does not require the existence of LTE whereas the 
second does and it is not clear how the results from the two methods can be compared. 

In 1-D systems with translation invariance, the 'conductivity', k, shows a power- law 
dependence on the size of the system; k ~ L a A question that has been investigated 
extensively is whether a is universal across a class of systems. The other related questions 
are about the number of such universality classes and the values of a characterizing them. 
Such an investigation has to take into account the difficulties posed by the considerations 
in the previous paragraph. Hydrodynamic mode-coupling theory predicts two universality 
classes corresponding to a = 1/2 and a = 1/3 for one-dimensional momentum conserving 
systems js, J]. Specific models of one-dimensional lattice systems with pair- wise nonlinear 
interactions between neighboring particles have been studied extensively using numerical 
and analytical methods [_2]. The analytical studies involve approximations at some stage, 
whereas numerical studies are handicapped by the system size and the time scales that can 
be achieved in computers which makes it difficult to give a definitive answer to the question 
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of the universality of a. Results of these studies agree with the hydro dynamic results that 
the generic value of a is 1/3; when the leading order nonlinearity is cubic a = 1/3 js|. 
However, for the so-called Fermi-Pasta-Ulam chain, where the nonlinearity is quartic, there 
is still no consensus; numerical simulations show different values of a 6|, 13] and analytic 
calculations, both classical and quantum, show a = 2/5 js-lQ]. 

Apart from the obvious theoretical interest in these problems there has also been a lot 



of experimental interest in studying heat transport at micro and nano scales 



11| . Recent 



experimental studies on nanowires such as Carbon nanotubes and Boron-nitride nanotubes 



show a clear breakdown of the Fourier law 
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15] . In order to attempt an understanding of 



realistic systems, it is essential to include all the modes of vibrations. Theoretical results, 
both numerical and analytical, mentioned above have been obtained for only the longitudinal 
vibrations of the chain. It is natural to investigate how the exponent a is affected by the 
transverse vibrations. 

Chains with both longitudinal and transverse vibrations have only been studied in the 



work of Wang and Li [16]. The Wang-Li model contains two kinds of infractions; a pair- 
wise interaction between the neighboring particles and another interaction that depends on 
the bending of the chain. The harmonic approximation of the model yields modes whose 
frequencies Q a {q) depend on the wave- vector q in the following manner. For small q, the 
longitudinal modes have £l x (q) cx q whereas the transverse modes £l y (q) oc q 2 . Wang and Li 
used the mode-coupling analysis (MCA) and molecular dynamics (MD) simulations to study 
this system. Their first important observation from MD simulations is that the phonon 
dispersions are strongly renormalized from their bare values, particularly the transverse 
phonon dispersion gets renormalized to become linear in q. For this reason they use an 
effective Hamiltonian method in which the form of the Hamiltonian and its parameters 
are taken from the MD simulations. For the conductivity exponent they obtain a generic 
value a = 1/3, but they also observe a value of a = 2/5 in their MD simulations which is 
attributed to crossover effects. When the parameter measuring the strength of transverse 
vibrations is sufficiently large, they also find that a = 1/2. 

In this paper, we study the quantum mechanical version of this model along the lines 
of flQ. In this procedure, one essentially ca.cn.ates the second order self- energy of the 
phonons self-consistently, which enables us to obtain the wave-vector dependent relaxation 
rate T(q). The thermal conductivity of the chain is then obtained through the use of Kubo 
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formula. The finite size of the chain is taken into account by putting on the time integral in 
the Kubo formula an upper cutoff proportional to the size of the chain jj]. This prescription 
yields a = 1 — 1/5, where 5 governs the small-q behavior of the relaxation rate, T(q) oc q s . For 
the Wang-Li model, we find that in the second-order perturbation calculation, the relaxation 
rate T y (q) for the transverse phonon diverges as q — > 0, whereas the relaxation rate T x (q) for 
the longitudinal phonons becomes constant. This makes the bare phonons invalid normal 
modes. These divergences basically arise due to the g 2 -dispersion of the transverse modes. 



This necessitates a self-consistent study of the phonon dispersion for this system [17H20]. 



Such studies have been done on several systems, including one-dimensional systems as test 



cases for the procedure 2l|-[23|. Here we follow the method of Gotze and Michel 18) which 



provides a nice connection between the Mori-Zwanzig mode coupling technique and the 
perturbation theory. We show that this procedure indeed renormalizes self-consistently the 
transverse phonon dispersion to be linear, in qualitative agreement with the numerical results 
of MD simulations [16]. We further note that the result on the renormalized frequencies is 
independent of the bare frequencies, but does depend on the fact that the cubic interaction 
vanishes linearly with q. Since the vanishing of the interaction with any of its wave-vector 
arguments is a consequence of the translation invariance, this result should be generic to 
translation invariant one-dimensional systems. 

We then perform a self-consistent analysis for the relaxation rates of the renormalized 
modes using the second order perturbation diagrams. This procedure leads to integral 
equations from which small-q dependence of the relaxation rates T a (q) can be extracted. The 
main results of this analysis are: T x (q) oc q 3 ^ 2 and T y (q) oc q 2 . Gratifyingly, this behavior 
matches the results derived by Wang and Li for the classical system using the mode-coupling 
analysis, which arrives at the results through a seemingly different mechanism. Finally we 
use these results for understanding the behavior of conductivity with the system length L. 

This paper is organized in the following manner. In Sec. II, we present the model 
and a calculation of the relaxation rates of the bare phonons. As mentioned above, this 
calculation shows that the modes are not well defined. Accordingly in Sec. Ill, we present a 
self-consistent calculation for the frequencies of the modes, which shows that both the modes 
should disperse linearly at small wave- vectors. In Sec. IV, we calculate the relaxation rates of 
the renormalized phonons using second order perturbation theory. In section V, we formulate 
equations for the relaxation rates using the self-consistent second order perturbation theory 
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and obtain their behavior at small wave-vectors. In Sec. VI, we turn to the calculation of 
thermal conductivity using the finite size procedure mentioned above. We close the paper 
with a summary and discussion of our results in Sec. VII. 



II. WANG AND LI'S MODEL 

This model consists of a chain of equal mass point particles with two kinds of interac- 
tions. First is a pairwise interaction between nearest-neighbor particles which depends only 
on the absolute value of the distance between the particles and the second is a three particle 
interaction between neighbors that takes into account the bending of the chain. There are 
two transverse directions, but due to the axial symmetry of the chain the two directions 
are entirely equivalent. Accordingly for understanding the qualitative role of the transverse 
vibrations on the relaxation of modes and transport, it suffices to consider just one trans- 
verse mode. In terms of the position r = (x,y) and momentum p = (p x ,p y ) vectors, the 
Hamiltonian is written as 

H = E + \ K r W + !) - r WI " «) 2 + K e com))} , (1) 

i 

where a is the lattice constant, cos(#(Z)) = — n(/ — 1) • n(Z) and n(Z) is the unit vector along 
r(l + 1) — r(l). The mass of the particles has been set to unity. This is a complicated 
many body problem. One assumes that a few lower-order terms in the Taylor expansion 
of the nonlinear potential is enough to give the generic physical picture. Also the low- 
temperature transport properties of the system are expected to be described by the low- 
energy modes of the system. So we set up the problem in terms of phonon modes. Denoting 
the deviation from equilibrium position of the particles by we have r(l) = (la + (p x (l), (j> y (l)) 
and p(Z) = (7r x (Z), 7T y (l)). Fourier transform along the length of the chain is defined by 

/(*) = -7^E e_a,rW /(0, (2) 



where N is the number of particles in the chain. The Hamiltonian, considering nonlinearities 



up to cubic terms only, can be written as [16] 

H = I E + h E ^(Mi, feaWU Wl(k) 



2 ^ 1 w tv ' 2! 

ki ki 



\= E ^(Mi, k 2 i 2 , k^lih^lih^h), (3) 
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where the summation indices in bold letters indicate that the summation is over all the 
polarization indices ii,i 2 ,-- which takes values from {x, y} and over all the wave-vectors 
ki,k 2 ,--, a notation we will follow hereafter. V^'s are symmetric under exchange of pairs 

invariance imposes two conditions on V; they vanish when any 
of the wave-vectors is set to zero and V n {kii\, .., k n i n ) oc A(ki + .. + k n ), where A(k) is 
unity when k is a reciprocal lattice vector and zero otherwise. Explicit forms for the Vs are 
obtained by expanding the potential term in Eq. (U) in Taylor series and taking the Fourier 
transforms. The harmonic terms are given by [16] 

V 2 {k 1 i 1 ,k 2 i2) = S^+kJi^Q^ki), 
Q x (k) = 2Q X \ sin(fca/2)|, Q y (k) = AQ y sm 2 (ka/2), 

n 2 x = K r , nl = K e /a 2 . (4) 

Note that the frequency of the longitudinal mode vanishes linearly with the wave-vector 
whereas for the transverse mode it vanishes quadratically. Henceforth we choose units in 
which a = 1. The reflection symmetry in the transverse direction permits only two kinds 
of cubic terms to be nonzero, these being V^{k\x, k 2 y, k 3 y) and V 3 (kix , k 2 x , k 3 x) . Both 
these interactions are proportional to sin(&i/2) sin(A;2/2) sin(&3/3). These forms are generic 
when the potentials are taken to be functions of absolute values of relative displacements 



(|r i+1 -r,|)(lfi| 



in [l0| i 



For the Fermi-Pasta-Ulam chain 10 in which only the longitudinal vibrations are con- 
sidered, the perturbation analysis is done around the bare phonon modes as they are well 
defined. But in this model we shall see that the bare phonon modes are not well defined, as 
the relaxation rates of the two kinds of phonons are not zero even for zero wave- vector mode. 
This requires that one works with the renormalized phonons as discussed below. We begin 
by calculating the relaxation rates of the bare phonons. For this purpose we first define the 
bare phonon operators 

in terms of which the quadratic part of the Hamiltonian becomes diagonal. The Hamiltonian 
can be written as 

H = ni(k)4(k)oi{k) + J2 Mhii, hi 2 , hi^U^U^UM, (6) 



3!ViV ki 
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FIG. 1: Diagrams for the second order contribution to the self-energy EJ- (q). The bare phonon 
Green's functions (k) for the longitudinal phonons and Dy (k) for the transverse phonons are 
represented respectively by the solid lines and the dashed lines. The diagrams (a) and (b) are two 
contributions for the longitudinal phonons, while the diagram (c) gives the sole contribution for 
the transverse phonons. 



where we have defined Ui(k) = o,i(k) + a\(—k) and V 3 (k\ii, k 2 i 2 , k 3 i 3 ) 



v 3 (kiii, k 2 i 2 , k 3 i 3 )/ y/ n h (ki)n i2 (k 2 )n i3 (k 3 ). 

For a system with two phonon polarizations, in general one defines the sin- 
gle particle Green's functions as the time-ordered expectation values Gij(q,r) = 



- T T 



ai(q,r)Oj(q,0) J, where the time dependence of operators is given by 0(r) = 
e rH Oe~ TH , t being the imaginary time. The angular brackets correspond to thermal aver- 
ages denned by (O) = Tr[e~ ,3H 0]/Tr[e~ l3H }. The G's are thus 2x2 matrices in polarization 
indices, but in the present case, the reflection symmetry in the transverse direction al- 
lows only the diagonal components to be nonzero. In the following we just write them 
as Gi(q,r). The perturbation series for Gj can be summed into the form Gj 1 (q,ico n ) = 
iuj n — £lj(q) — Tij(q,icu n ), where u^s denote the usual Matsubara frequencies and Hj(q,iu n ) 
is the self-energy. The relaxation rate is defined as Vj(q,u) = —Qz~2j(q,iu n — >■ u + i0 + ), 
and the on-shell relaxation rate is denoted by Tj(q) = Tj(q,flj(q)). We first evaluate the 
relaxation rate up to second-order in the perturbation series. The second-order diagrams for 
S x that contribute to the relaxation rate are given in Figs. [T] (a) and (b). Their contribution 
is 

£i 2 W) = ~ E \%{-^Mhk 2 i)\ 2 Dt\k U T)Df\k 2 ,T), (7) 

ki,k2,i 
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where Df*\q, r) is the following unperturbed Green's function. 

Df ) \q,T)=Gf\q,r) + Gf\-q,-r) = - £ e - s ^n s (^(<?)), (8) 

se{+,-} 

where n_(f2) = n(f2) is the Bose factor and n + (Q) = 1 + n(Q). The relaxation rate is given 
by 

4 2) (?) = ^(l-e-^ (9) ) Yl \%{-wMiMi)?n s M{h))n s M{h)) 

ki,k2,si,S2,i 

x5{Q x {q) - stttiih) - 8 2 Qi{k 2 )). (9) 

The sum is simplified by two conditions: the A-function contained in the potential V3 and 
the ^-function. Therefore, we can write 

x5(fi x (g) - siOi(fc) - s 2 fii(9 - k)) 
= i(l-e"^W) £ |V 3 (- 9 x,^g-^)P g _ 



51,S2 



where £;* = fc*(g, si, s 2 , i) is the solution of fc(q, k, si, s 2 ) := fia:(<z) — Sifij(/c) — s 2 f2j(g — fc) = 
and the Jacobian Ji(q, k* , S\, s 2 ) = (d/dk)fi(q,k,si,s 2 )\k=k*- If the sum on the RHS of 
Eq. (fTOj) is finite we expect T x 2 \q — > 0) oc ^(9), as |t^(— gx, ki,q — ki)\ 2 oc We first 

note that the term corresponding to i = x (Fig. 1(a)) which involves only the longitudinal 
modes gives zero contribution as the condition / = can have only trivial solutions at which 
the self-energy itself becomes zero [5] . 

The second term corresponding to i = y (Fig. 1(b)) needs to be examined for singular 
behavior. The source for singularities are the Jacobian J and the number factors which for 
low frequencies have the form n(Q — > 0) ~ (/3f2) -1 . For this term we have transverse mode 
frequencies which vanish quadratically with the wave- vector and hence are more singular. 
Considering s\ = + = s 2 term, the small-g solution is obtained to be k*(q) ~ y/\q\ and the 
Jacobian J(q, k*) ~ y\q~\- This yields a contribution to T x (q) which diverges for small q at 
finite temperatures and we have the leading behavior T x (q — > 0) ~ |g| -1 ^ 2 - 

The second-order diagram for H y is given in Fig. 1(c). The yields the relaxation rate as 

T y (q) = I(l-e-^W) JdkY,\V(-qy,kx,q-ky)\ 2 n Sl (n x (k))n S2 (n y (q-k)) 

sis 2 

x6(tt y (q) - SiQ x (k) - s 2 Q y (q - k)). (11) 



Proceeding as in the analysis of Eq. ( flUj) . we see that the term corresponding to s\ = +,s 2 = 
— is finite for q — > 0. 



III. RENORMALIZATION OF PHONONS 

The above analysis shows that the bare phonon modes are not stable under scattering 
and thus are not good zeroth order approximations in the perturbation theory. Therefore 



we follow the mode-coupling method as formulated by Gotze and Michel 18] to study how 



the bare phonon dispersions get affected by the nonlinearities. For completeness we first 



present a brief outline of the mode-coupling method [18J. The procedure requires evaluation 
of scalar products defined in the space of the dynamical operators (denoted by Oi, 2 , ■■) as 
follows. 

(0 1; 2 ) = -<Ol;0 2 »*=o, 
« o i; 2 » z = ± J dtQ(±ty zt ([O^t), O 2 (0)]> , Qz % 0, (12) 

where 0(t) = e lHt O = e Oe~ lHt . For us the key quantity is the static suscepti- 
bility (q) = ((f>i(q), 4>j(q)). It is easily seen that in the harmonic approximation 
Dijiq) = &i(q)fiij- The renormalized frequencies are the eigenvalues of the matrix D(g) 
for the nonlinear problem. It follows from the definition and the space reflection symme- 
try of the Hamiltonian that D(g) is a real symmetric matrix which is an even function of 
q. Therefore, D is diagonalizable with a complete and orthonormal set of eigenvectors e a 
obeying 

Y J D ij {q)^{q)=ujl{q)e^{q), (13) 
3 

where the eigenvalues obey oo a (q) = u a (—q) due to time-reversal symmetry. 

In order to develop equations for these scalar products, we utilize the identity, 
(O u n0 2 ) = ([Ol,0 2 ]) which follows from Eq. (TJJ). The LHS of this equation 
can be calculated using HO = [H,0]. For further calculations we define Ai(q) = 
l^k i V(qh-kiii,-k2i2)<Pi 1 (ki)<Pi 2 (k 2 )- Taking (i) O x = <p(q),0 2 = ir(q) and (ii) O x = 
A(q),0 2 = vr(g) respectively, we get the following equations which are exact. 

DiM) = V(~qi, qj) + ^ (A(q), <j> k (q)) D kj (q), (14) 

k 

J2(Mq), MqW(-qk, q i) = -(A( q ), am))- (ib) 



We write these equations in the matrix form, suppressing the g-dependence, as 



D = V[V+ (A, A)PV 



(16) 



where (A,A) y (g) = (Ai(q),Aj(q)) = - < Ai(q)\Aj(q) > 2=0 and V y (g) = V 2 (-gi,gj)- To 
get a closure for this set of equations, we need to make some approximations. We do this in 
a self-consistent manner by expressing the higher order correlation functions (Ai(q), Aj(q)) 
in terms of D matrices. The expression for (A, A) is is given by 



< Ai(q); Aj(q) » 2 = ^ V 3 (-qi, -Mi, -M2) V 3 (gj, -k 3 i 3 , ~hu) 



ki 



(17) 



The expectation contained in the above equation is evaluated using an effective quadratic 
Hamiltonian given by 

#0 = ^^(^w + ^E^w^w^w. 



Ai 



A,tt 



w a (/c)6l (fc)&a(*0 + - 



(18) 



where the operators 6 a (A;) are the bosonic operators defined in the usual manner as 
b a (k) = Y,e«{k) \y/2u) a (k)(f)j(k) + i\/2/uj a (k)7ij(k) 



(19) 



Here the frequencies oo a (k) and the polarization vector ef(k) are defined in Eq. f|T3|) . and 
are to be determined self-consistently. 

Under this approximation using Eq. ( [T2|) we get 

(A(a) A (a))-2^G..(a 1 J l + »(!) +n(2) n(l)-n(2) 
(^g),^(g))-2^(g,l,2) j w(1) + w(2) + w(2) _ w(1) 

^Ml^li^^ (20 ) 

ki 

where we have used the compact notation 1 = (ki,ai). Further w(l) = w ai (fci) and n(l) = 
n(u;(l)). Eqs. ffTB"]) and fl20l together provide us the self-consistent equations for determining 
the renormalized frequencies u a {q)- In order to obtain their solution, we first iterate them. 
Starting with a diagonal Dij(q) = 8ijdf(q) in the RHS of Eq.( l20]) . we have e" = 5 ia and 



we 
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see that only the diagonal terms survive. So this procedure gives a diagonal D back if our 
starting D is diagonal. Putting this information in Eq. f|T6|) reduces it to 

D (a) = ^ 

(A A) (q) - - Ij^Cg^fciVM)]^ f 1 + n(u)j(ki)) + n(u)i{k 2 )) | njujjjki)) - n{uj(k 2 )) \ 

, a as / \ 1 |V3(gy, fc 2l/)| 2 / 1 + n{u x {k 1 )) + n(u y (k 2 )) n(w s (ki)) - n{u y {k 2 )) 



2^ u) x (ki)ujy(k 2 ) { oa x (ki) + u y (k 2 ) u y (k 2 ) - u x (h 



In order to analyze the low wave- vector (q — > 0) behavior of D(q), we simplify the equations 
by noting that 

V 3 (qi, 1, 2) oc A(g + k\ + fc 2 ) sin(g) sin(A;i) sin(/c 2 ) > A(g + k x + fc 2 )g sin 2 (/c 1 ), 
ra(u; a (fci)) - n(u a (-kl - q)) q ^o d 

7 — T. \ 7T\ y ~~E~ n ( U} <x{ki)). (22) 

w a (- k\ — q) — UJ a {ki) OLO 

We first consider (A, A) xx (q), whose expression contains two terms corresponding to scat- 
tering by longitudinal (i = x) and transverse (i = y) modes respectively. We may first 
suppose that the renormalized frequencies have the same g-dependence as the bare fre- 
quencies, oj(q) oc 0(5). Then the use of Eq. f l22|) in Eq. f l2T|) shows that, as q — > 0, the 
longitudinal contribution to (A, A) xx (q) behaves as q 2 whereas the transverse part gives rise 
to a divergence basically due to g 2 -dependence of the transverse frequencies in the denomi- 
nator. Similarly the expression for (A, A) yy (q) also involves transverse mode frequencies in 
the denominator, which again give rise to a singular dependence on q. 

If on the other hand we assume u a (q — > 0) oc q, all the integrals converge and we get 
D a a(q — > 0) oc q 2 . Thus we have a self-consistent solution for D aa (q) in the q — > limit, 
which is also diagonal. This is exactly what is seen in the MD simulations of the classical 
system with the Wang-Li interaction [if]]. From this analysis it is clear that the linear 
solution obtained here is independent of the form of bare frequencies. Thus the linear in- 
dependence of the frequencies in one dimension seems to be a generic result. For further 
analysis we take a specific form for renormalized phonon frequencies, namely, 

Ui(q) = UiU)(q), u(q) = 2| sin(g/2)|, u y = r]u x , (23) 

which agrees with the self-consistent solution in the q — >■ limit. 



11 



IV. RELAXATION RATE OF RENORMALIZED PHONONS 



We now study an effective Hamiltonian given by Eq. fl6]) and Eq. (J5J) in which fl(q) is 
replaced by u(q). The second-order relaxation rate for the renormalized modes are as before, 

4%) = i(l- e -^(9)) J dk^\V 3 (-qx,ki,q-ki)\ 2 n sl (ui(k))n S2 (uJi(q-k)) 

IS1S2 

xS(u x (q) - SiUi(k) - s 2 uji(q - k)) (24) 
Tf\q) = l(l- e -^)) / dkJ2\V3(-qy,kx,q - ky)\ 2 n Sl (oo x (k))n S2 (oo y (q - k)) 

x5(u y (q) - siuj x (k) - s 2 u y (q - k)), (25) 
V 3 (hh, hh, k 3 i 3 ) = V(kii u k 2 i 2 , k 3 i 3 )/ y/ u h (h)uj i2 (k 2 )u i3 (h) (26) 

(2) 

First we estimate the low-g form of T y (q). This proceeds exactly as the analysis of Eq. (imp . 
We have 

rf( 9 ) = i(i-e-""«<") Y, ^'uu'vl ^'"' vMyiKMi - *•))■ (27) 

where fc* are the solutions of / = ^(g) — Siu^fe) — s 2 u y (q — k) = and the Jacobian is 
given by J = df/dk\k=k*- Since both uj Xjy (q — )■ 0) ~ g the factors of g in the denominator 
that arise due to number factors are canceled by the \V\ 2 term. Therefore if the Jacobian 
J(k*) does not go to zero at some k*, we get T y 2 \q — > 0) oc oo 2 (q). 

There are four cases for the energy-momentum conservation conditions due to four 
possible values of s± and s 2 . These have to be analyzed separately. We choose the 
Brillouin zone [0, 2ir]. For si = + = s 2 , the energy conservation condition becomes 
oo(q) — u(k)/i] — u(q — k) = 0. When r\ — 1 there are two solutions k = and k = q. 
But these correspond to no scattering and have to be excluded from the perturbation terms. 
Thus as noted earlier the relaxation rate within this approximation is zero js|. When r\ > 1 
nontrivial solutions are possible. Fig. [2(a) shows plots of uj(q) and u(k)/i] + u(q — k) with 
respect to q for a representative set of values for k and rj. This figure shows that, when 
77 > 1, k is a solution for the energy conservation condition oo(q) = u(k)/r] + u(q — k) when 
the value of q corresponds to one of the intersections of the two graphs. For small-g we 
get a solution k ~ q with a finite value for the Jacobian of the argument of the 5-function. 
Therefore we see that the contribution to T y 2 ^ (q — > 0) from this case is proportional to q 2 if 
i] > 1 and zero otherwise. 
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FIG. 2: Figures demonstrate how the energy conservation condition in Eq. (|25j) can be satisfied. 
The figure on the left shows ui(q) and uj(k)/r) + u)(q — k) in the same graph for k = 0.5, r/ = 1.5 
and the figure on the right shows u:(q) and u(k)/r] — oo{q — k) for k = 0.5, rj = 0.5. Dotted vertical 
lines mark the values of k which is the solution for q values corresponding to the intersections of 
the two graphs. 

For si = +, S2 = — , the 5-function condition is u(q) = u(k)/r] — co(q — k). Fig. [2(b) 
shows LHS and RHS of this condition on the same graph. We see that for rj < 1 we have 
a nontrivial solution and in the small-g limit we find a solution k ~ q where the Jacobian 
also is finite. Therefore, for this case, the contribution to T^y \q — > 0) is proportional to q 2 
for rj < 1 and zero otherwise. 

Next we take Sj = — , S2 = +■ The energy conservation is given by oo(q) = u(q — k) — 
u(k)/r]. We have a nontrivial solution for r] > 1 which for q — > is k ~ q. Again we find the 
contribution to Ty \q — > 0) is like q 2 for r] > 1 and zero otherwise. The fourth case Si = S2 
has only trivial solutions. Therefore, we have 

r [ y\q^0)~ q 2 if 77 ^ 1. (28) 

(2) 

Now we consider the relaxation rate T x . It has contributions coming from both trans- 
verse and longitudinal modes corresponding to i = y and i = x in Eq.( [2^1) . The i = y 
term again involves interaction between two ?/-phonons and one x-phonon as in the case 
for and its low wave-vetor behavior can be calculated in a similar way. As in the cal- 
culation of 1^ , for rj ^ 1, there are nontrivial solutions satisfying the ^-function except 
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for s\ = S2 = — 1. Thus the corresponding contribution to T x (q — > 0) is proportional to 
q 2 . The i = x term involves three x-phonons and the contribution to the relaxation rate is 
zero because the 5-function condition sin(g/2) — Sisin(A;/2) — s 2 sin((g — k)/2) = has no 
nontrivial solution. Therefore 

if^O)^ 2 ifr ^ L (29) 

Thus, the renormalized phonon modes at low wave-vectors are slow decaying and well- 
defined to perform the perturbation analysis. 



V. SELF-CONSISTENT RELAXATION RATE 



The phonon relaxation rates of a chain with cubic nonlinearity but with only the longi- 
tudinal vibrations was studied earlier by us [5|. For this mentioned in the above 
paragraph, the relaxation rate is zero in the second order as the two conservation conditions 
for energy and wave- vector cannot be satisfied simultaneously. To overcome this difficulty we 
performed a self-consistent approximation which consists in replacing bare phonon Green's 
functions in the second-order diagram (FigffJ^a)) by full Green's functions. This relaxes the 
strict energy conservation condition and allows for nonzero solutions for the relaxation rate 
through a solution of an integral equation. Our aim in this section is to perform a similar 
analysis for the present case. For this purpose, one uses the spectral representation of the 
Green's function, 

A(g,r) = - Yl / ^#,e)n s (e)e- SfT , (30) 
with the spectral function is given by 

m ' e) = mJItim' (31) 

where Ui(k,e) = Ui(k) — Sf (g, e) and Sf (g, e) = 3ftSj(g, iu n — > e + i0 + ). If we now replace 
D^(q,r) by Di(q,r) in Eq. ([7]), the corresponding equations for relaxation rate, Eqs. ((2] 
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and (125]) . are replaced by a set of self-consistent equations for T x and Y y . We get 




k.s.i 



xB i (s 2 k 2 ,s 2 U) x (q) - s 1 s 2 e)n 8a (s 2 fjJ x (q) - sis 2 e), 



r tf (?) 




k.s 



xB y (s 2 k 2 , s 2 u y (q) - SiS 2 e)n S2 (s 2 o;j / (g) - Sis 2 e). 



(32) 



We follow Ref [5j in the analysis of these equations. We note that for small q, the relax- 
ation rate Ti(q) is very small and the spectral function Bi(q,e) can be approximated by 
a narrow peak at e = u)i(q,Ui(q)) =: u)j(g). Also, d)j(g — > 0) oc Further we may 

replace B i (p,e)f(e) ~ B i {j>,e)f{u i {p)) where /(e) is a smooth function of e. We use this 
approximation for the number factors n(e) in Eq. (1321 . We also replace a)j(g, e) ~ u>i(q) and 
Tj(g, e) ~ Tj(g) in the expression for B^q, e) to be used in Eq. (|32|) . After these approxima- 
tions, only the spectral functions depend on e in Eq. fl32l) and the integral can be performed 
using 



We note that tJi(q) = Ui(—q) because it depends on q through cu(q) which is even in q and 
through V 3 factors which enters the perturbation terms as even function of q only. Then, 
Ti(q) = Ti(—q) and n(q) = n(u(q)) = n(—q). We can write T x (q) = T x (q) + T x (q), with 




(T 1 + r 2 )[(a-b) 2 + (T 1 -T 2 ) 2 ] 
{{a-by + Tl + TlY-ATlTl ' 



(33) 



1 (1 _ e - M) J^nMk^nMq + k))\V 3 (-qx, —ki, -q -ki))\ 2 



xS ++ (qx, ki,q + k i) + e^ xi - q) n(uji(q - k))\V 3 (-qx, ki,q-ki)\ 2 



1%) 




(34) 



xS + -(qx, ki,q + ki), 



(35) 



where 



Sssssifah, k 2 i 2 , k 3 i 3 ) 



(r i2 (fc 2 ) + r i3 (fc 3 ))[xL 2 ( fc A>fc2»2,M3) + (T t2 (k 2 )-r l3 (k 3 )) 2 ] 

[X 2 sis2 (kii h k 2 i 2 , k 3 i 3 ) + Tl(k 2 ) + Tl(k 3 )} 2 - 4T? (* a )I* 



.(36) 
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Here we have introduced the definition: Xsis 2 (kiii, k 2 i2, ^3*3) = ^i 1 {ki) + S2^i 2 {k2) + s 3 uii z {k 3 ). 
Similarly, for the transverse modes we get 

T y (q) = 1(1 - e-^W) ^ {nK(A;))nK(g + fc))|V- 3 (-g2/, -fcr, g + ky)\ 2 S ++ (qy, kx, q + ky) 

k 

+e^ q) n{u x {k))n(ujy{q - k))\V 3 (-qy, kx, q -ky)\ 2 S—(qy, kx, q-ky) 

+n(u x (k))[l + n(u y (q + k))]\V 3 (-qy, -kx, q + ky)\ 2 S + ^(qy, kx, q + ky) 

+ [1 + n(u x (k))]n(u y (q - k))\V 3 (-qy, kx, q - ky)\ 2 S-+(qy, kx, q - ky)}. (37) 

Equations (I3~4j) . (1331) and (13T1) have to solved self-consistently to obtain the small wave- 
vector behavior of the relaxation rates. We show that these equations have a solution 

T x (q ->• 0) ~ g 3/2 , T y (q -> 0) ~ g 2 . (38) 

The key feature that enables the solution of these equations is that for small q, each of the 
equations can be written in the following schematic form 

r,(g 0) oc /H 2 (<?) y dk F(k, q). (39) 

Here one factor of Ui(q) comes from the prefactor, which is basically a detailed balance factor 
and the other from the factor of | V3 1 2 . The integrand F(k, q) is not known as it as it depends 
on rj(fc)'s. Our strategy is to examine these integrands in Eqs. (134|) . (1331) . (13T1) at q — 0, by 
first assuming that Ti(k) oc k 2 . If the integral turns out to be finite, clearly Ti(q) oc q 2 is 
a self-consistent solution. If on the other hand the integral diverges, it presumably behaves 
for small q as q~ x . So we assume that T(k) ~ k s with 5 < 2, and examine the integral for 
small q. A solution is achieved if a self-consistent value for 5 can be obtained. In examining 
the singularities of the integrand in Eq. (1391) . a further simplification comes from noting 
that the Bose factors n(k) give rise to factors l/uji{k) for small frequencies, but these are 
canceled by the similar factors in interaction term \V\ 2 . So all the singular behavior of the 
integrands come from the S'-factors introduced in Eq. (13~6l) . 

First, let us consider Eq. (13T1) . the equation for T y (q). At q = 0, the S'-factors could 
give rise to singularities at k = if we set T(k) cx k 2 . For 77 7^ 1, in all the four terms 
in Eq. [371 x{^)kx,ky) ~ k for small k; for example, in the third term, x + _(0, kx, ky) = 
2(1 - n)wn(k/2). Using Eq. ([38]) one finds that S S1S2 (0, kx, ky) ~ P*; 2 /A; 4 = for small 
k. Since the integral is finite at q = 0, r y (g — >■ 0) ~ g 2 . We next consider T^g) which has 
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been split into two terms r^(g) and T^(q) in Eqs. (|34|) and (1351) . The analysis of r^(g) 
proceeds in a similar way as x(0, ki) ~ /c as — > in both the terms of Eq. fl34|) and we 
get r£(q) 0) ~ g 2 . 

Something different happens in Eq. 035]) for (g)- It has again two terms for i — x,y. The 
term 2 = x involves only the longitudinal mode, and this situation is identical to the FPU-a 
problem discussed in Ref. jsj. We note that x+-(<2S kx, q + kx) ~ q[l — cos(A;/2)] oc qk 2 as 
q, k — > 0. Using Eq. f l38j) . for small g and k, we get 5* + _(gx, /cx, q+kx) oc [r x (g)+T :r (g+A;)] _1 . 
At g = 0, the corresponding integral is divergent. But the small q divergence of the integral 
can be obtained by noting that the major contribution to the integral comes from small 



k and smal 



following 



5] 



q + k region. Denoting the corresponding contribution to by r 1; we get, 



Jo T 1 (q + k) +T 1 (k) 



As is easily checked, Eq. ( 140]) is solved by taking Ti(k) cx k 3 ^ 2 . For the i = y term 
X+-{qx,ky,q + ky) ~ g[l — rjcos(k/2)] for small g and k. The leading order behavior 
of £ in this limit is obtained using T x (q — > 0) ~ g 2 . We see that S + ^(qx, ky, ky) is not 
singular at k — for finite g, S + _(qx, Oy, qy) ~ g 2 g 2 /g 4 ~ g°. Thus the limit g — )■ of this 
sequence is finite and the corresponding contribution to (g — >■ 0) is proportional to g 2 . 
We can now go back to the earlier equations and check if by setting T x {k) oc fc 3 / 2 . the earlier 
results for T y (q) and T^(g) are affected. The examination of the integrals shows that they 
remain finite. Therefore to the leading order, T x (q — > 0) ~ g 3 / 2 + 0(g 2 ) and T y (q — > 0) ~ g 2 
establishing the solution given by Eq. fl38|) . 



VI. THERMAL CONDUCTIVITY 

In this Section, we calculate the thermal conductivity of the system using the Kubo-Green 
formula. Defining the thermal Green's function for the total current operator J as, 

W(r) = ~<r T [J(r)J(0)]>, (41) 



the Kubo-Green formula for thermal conductivity can be written as 



24] 



« = -ilim lim ^(iu^u + iS) 

T w->-0 7V->oo OJ 
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where oo n are the Matsubara frequencies. The current operator, up to leading order, is given 
by 

J = -J2Uk)bt(k)b a (k) (43) 

k.a 

where j a {k) = v a (k)fl a (k), v a {k) = (d/dk)u a (k). Further, we neglect the cross- correlation 
between the x- and y- components of the current. Under the approximation used in [h]] that 
vertex corrections to the two particle Green's function Eq. ( 1411 has the same wave-vector 
dependence as j(k) for small k, one arrives at the kinetic theory formula for conductivity, 



K 



£ / ^l(q)T a (q)C a ( q ), (44) 



where r a (q) = r~ : (g) is the relaxation time of single particle Green's function and C a (q) = 
i(— ■^n(u)) ul=ulct ^ q )u'^(q) is proportional to the mode specific heat. We note that C(q) is 
well behaved for small q and any possible divergence can come from r(q) only. 

The formula given by Eq. f|44l is for an infinite chain and this quantity is divergent 
because r a (q) has a non integrable singularity at q = 0. To obtain the conductivity k(N) for 
large but finite size chain of N particles, we follow the prescription used by other workers in 
the area [2|. In this procedure one rewrites Eq. (1441) as a time integral in the following way. 

K = J2 [°°dtK a (t), K a (t) = [ %l{q)C a {qy~ Ta{q)t (45) 
a Jo J 271 

Since T a (q — > 0) ~ q Sa , we get K a (t — > oo) ~ t~ l l 5a . One argues that when the boundaries 
of the system are connected to thermal baths, the correlation K(t) dies after a time scale 
of order N due to interaction with the baths, i.e. the typical time it takes for the modes to 
travel across the chain. Using in the time integral in Eq. (145]) such a cutoff yields 

k(N) ~ O^N 1 - 1 ' 8 *) + 0(N l - l ' & y) ~ 0(iV 1/3 ) + 0(N 1 ' 2 ) ~ N 1 / 2 . (46) 

It is to be noticed that the leading order divergence arise from the term corresponding to the 
2/-modes. If the conduction was normal, conductivity would be finite and the contribution 
from these modes would be negligible compared to the contribution from the x-modes when 
the strength of the interaction giving rise to transverse oscillations, Kq, is small. But here 
the conductivity itself is a divergent quantity and its divergence is determined by the slower 
relaxation of the transverse modes compared to the longitudinal modes. 
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VII. SUMMARY 



We have studied thermal transport in a chain of coupled particles in which particles can 
have longitudinal as well as transverse vibrations. Particles interact via anharmonic poten- 
tials proposed by Wang and Li 16|. For this model, the mode frequencies in the harmonic 
approximation have the following dispersions at small wave-vectors: for the longitudinal 
mode Q x (q) oc q; for the transverse mode Q y (q) oc q 2 , where q denotes the wave- vector. 
We first show that due to quadratic dispersion of the transverse mode, the relaxation rates 
for both the modes become too large to leave them as meaningful quasi-particles. So we 
calculate their renormalized frequencies taking into account the anharmonic interactions in 
a self-consistent approximation. We find that this procedure yields the linear dispersion for 
both kinds of modes irrespective of their bare forms. We then perform a calculation of the 
relaxation rates for the renormalized phonon in a self-consistent approximation based on 
the second-order perturbation diagrams. This calculation yields that for the longitudinal 
phonons, the relaxation rate T x (q) oc q 3 ^ 2 , while for the transverse mode T y (q) oc q 2 . We 
then use these results to calculate the thermal conductivity k(N) of a chain of N particles. 
The finite size is incorporated by using a cutoff of the order of N in the time integral over 
current- current correlation function that occurs in the Kubo formula for the conductivity. 
This yields a contribution to k(N) from longitudinal phonons to be 0(N 1 ^ 3 ), that from 
transverse phonons to be 0(N 1 ^ 2 ). The latter contribution dominates giving k(N) oc N 1 / 2 
for large N. 
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